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Abstract 

Most finite element methods for solving time-harmonic wave - propaga¬ 
tion problems lead to a linear system with a non-normal coefficient matrix. 
The non-normality is due to boundary conditions and losses. One way to 
solve these systems is to use a preconditioned iterative method. Detailed 
mathematical analysis of the convergence properties of these methods is 
important for developing new and understanding old preconditioners. Due 
to non-normality, there is currently very little existing literature in this 
direction. In this paper, we study the convergence of GMRES for such 
systems by deriving inclusion and exclusion regions for the pseudospec¬ 
trum of the coefficient matrix. All analysis is done a priori by relating the 
properties of the weak problem to the coefficient matrix. The inclusion is 
derived from the stability properties of the problem and the exclusion is 
established via field of values and boundedness of the weak form. The de¬ 
rived tools are applied to estimate the pseudospectrum of time-harmonic 
Helmholtz equation with first-order absorbing boundary conditions, with 
and without a shifted-Laplace preconditioner. 


1 Introduction 

Several different strategies for discretizing time-harmonic wave propagation prob¬ 
lems using finite elements have been proposed in the literature. For typical 
problems, most of these strategies lead to a linear system with a large, sparse, 
indefinite and non-normal coefficient matrix. The indefiniteness is due to the 
wave-nature of the problem and the non-normality arises either from losses or 
truncation of infinite domains to finite ones. The large size of the system is due 
to the number of degrees of freedom required to resolve an oscillating solution. 
Because of their properties, the linear systems related to time-harmonic wave 
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propagation problems are difficult to solve. Memory is an issue with direct 
solvers and lack of efficient preconditioners with iterative ones. 

In order to develop new preconditioners and to understand old ones, it is 
important to know their effect on the convergence properties of the applied 
iterative method. Unfortunately, the convergence of iterative methods for lin¬ 
ear systems with a non-normal coefficient matrix is a difficult subject of study. 
When the non-normality is significant, the iterative properties can be very dif¬ 
ferent from the ones indicated by eigenvalues, EDiini. Similar difficulties are 
met with other properties related to the non-normal matrices, e.g., behavior 
of matrix exponentials cannot be predicted by eigenvalues, m- Determining 
when the non-normality has a significant impact to iterative properties is com¬ 
plicated. The simplest way to estimate the impact is to compute one of the 
commonly used scalar measures of non-normality, e.g., — 74*A||||A||“^, 

the conditioning of eigenvectors or the conditioning of individual eigenvalues, 
[24] . However, except for the first one, these measures are not computable for 
large matrices. In addition, they can vary considerably even for relatively small 
systems m- 

The convergence of preconditioned iterative methods has been extensively 
studied in the context of finite element methods, [22]. However, majority of the 
research has focused on real valued symmetric positive definite problems. The 
finite element discretization of these problems also leads to symmetric positive 
definite linear systems, which are solved using the preconditioned conjugate 
gradient method (PCG). The aim in the analysis of these methods is to estimate 
the convergence rate before computations. Only few of the existing works deal 
with indefinite linear systems, iHiiaiaisiiis], and even fewer with non-normal 
indefinite ones, [ssiiiaES]. 

Most preconditioners for finite element discretizations of elliptic weak prob¬ 
lems have been analyzed by using the abstract framework of Schwarz methods, 
[22] . This framework is based on studying the properties of the underlying weak 
problem instead of the linear system. The convergence of PCG is related to the 
weak form via Rayleigh quotients. Such analysis is done in the inner product 
induced by the bilinear form. As Rayleigh quotients are the first step in the ex¬ 
isting analysis, it does not carry over to complex valued, indefinite, non-normal 
linear systems. Such systems require a different set of analytics tools. 

There currently exists three different ways to analyze iterative properties of 
a non-normal matrix [5]: to study the field of values (FOV), pseudospectrum, 
or to include conditioning of eigenvectors to the convergence estimates. For 
time-harmonic wave-equations, estimating eigenvector conditioning before the 
matrices are constructed seems to be complicated and thus this approach is not 
suitable for our purposes. FOV has been applied to analyze the preconditioned 
time-harmonic Helmholtz equation e.g. in [12]. However, FOV is always a 
convex set containing all eigenvalues of the matrix. As we will see, the spectrum 
of the problems we are interested in curls around the origin making FOV based 
methods unsuitable for our purposes. In contrast, the pseudospectrum can be a 
non-convex set and as we will show it can be estimated a priori, making it the 
best option of the three for this work. 


2 


In this paper, we study pseudospectrum as a tool for relating the properties 
of the weak problem to the convergence of the GMRES method. We derive 
convergence estimates for GMRES by establishing inclusion and exclusion re¬ 
gions for the pseudospectrum. The exclusion region is derived from the stability 
estimates of the weak problem and the inclusion region is based on then relation 
between pseudospectrum and EOV. In several cases, an inclusion for EOV can 
be easily obtained based on continuity properties of the weak form. All analysis 
is done a priori, so that the regions can be obtained without constructing the 
actual matrices or performing computations with them. The derived bounds are 
explicit in the relevant parameters of the problem, e.g., mesh size, wave-number 
and the losses. The presented analysis relies on general properties of the weak 
problem, stability and continuity so it is possible that it can be applied to other 
preconditioners and problems. 

The paper is organized as follows. We begin with some preliminaries and 
proceed to give estimates relating pseudospectrum to a weak problem. After 
establishing these abstract results, we apply them to three example problems. 
We begin the examples by considering the Poisson problem, which is included 
for easy reference on what kind of information the derived estimates can de¬ 
liver. Then we apply the presented tools to time-harmonic Helmholtz equation 
with and without a shifted-Laplace preconditioner. We end the paper with a 
discussion of the presented material. 

2 Preliminaries 

Our model problem is: Eind u G V such that 

a{u,v) = L{v) V'i; G E, (1) 

where V is some finite element space, a(-,-) : V x V ^Cisa sesquilinear 
form and !/(•): E ^ C an antilinear functional. The finite element space E is 
spanned by basis so that every function u E V admits the representation 

N 

u = 


in which the vector of coefficients Xu G C^. Problem m leads to the linear 
system 

Axu = h, 

where A G G C^, Aij := a{ipj^ipi) and hi := L{(pi). Hence, the 

sesquilinear form and the matrix A are related as 

a{u,v) := xlAxu, (2) 

where * - is the conjugate transpose. The properties of the matrix A will depend 
on the properties of the sesquilinear form and the basis functions via the above 
equation. 


3 


We will describe the actual problem and discretization in detail in Section 
4. For now, let us note that when the sesquilinear form a is related to the time- 
harmonic Helmholtz equation with absorbing boundary conditions, the matrix 
A can be very large. This is due to the facet that the finite element mesh size 
has to be sufficiently fine before finite element method can produce accurate 
results, see Ham]. Typical engineering rule of thumb is to use ten degrees of 
freedom per one wave-length. For example, a cube for which each dimension 
is ten wave-lengths long requires one to use 10 ^ degrees of freedom, this is, 

= 10 ^ or larger. 

In the following, we assume that problem m has a unique solution and ad¬ 
mits some kind of a stability estimate. Stability estimates are typically derived 
under additional assumptions on the domain and the antilinear functional L. 
In general, the functional L can be from the space V' = {f : V C \ f gH*}, 
where V* is the dual space of V. As such functionals can be quite badly behav¬ 
ing, stability estimates are often derived under the assumption L G IF', where 
V C IF. In this spirit, we make the following assumption. 

Assumption 2.1. Let W be a Hilbert spaee, V C IF, L e W' and u e V 
be the unique solution to problem m- Then there exists a eonstant Cs > 0 
independent of u and L sueh that 

||w|| < Cs\\L\\w' (3) 

where || • || is a norm on V and || • ||w' •= sup{ \L{w)\ \ w eW and ||u;||w = 1 }• 

The pseudospectrum of a matrix A G Ae(A), is a family of sets 

depending on a parameter e > 0. The sets in the family are defined as 

V(^) := { 2 e C I \{zl - A)~^\>e~^} , 

in which | • | is the standard spectral norm. When the matrix (zI — A) is singular, 
we define \{zl — A)~^\ = oo. The notation | • | is also used for the Euclidian 
norm of a vector. Clearly, the pseudospectrum can also be characterized as 

Ae(A) := { z eC \ (TminizI - A)< e} , 

in which we denote the smallest singular value of a matrix B G as 

^min (-^) • 

The pseudospectrum was independently proposed by several authors as an 
extension of the spectrum, suitable to study the properties non-normal matrices, 
m- The pseudospectrum has been extensively studied in the literature, see e.g. 
niiiiin]. In the following, we write A^{A) = A^, when the matrix A is clear 
from the context. 

In the derivation of the inclusion region, we take advantage on the relation 
between FOV and pseudospectrum. The FOV of a matrix A G is defined 

as the set 

FOV {A) := I ^^ I X G C^and x 7 ^ 0 | . (4) 
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The set FOV{A) is convex, compact and contains all eigenvalues of A. As 
we will see, coarse inclusion for FOV{A) can be obtained by using it’s close 
relation with the sesquilinear form. We postpone stating the relation between 
pseudospectrum and FOV to Section 3, where we have introduced sufficient 
notation for proving it. 

Both pseudospectrum and FOV can be related to convergence of the GMRES 
method, 13 uni HI]. The approximation error for the solution Xi generated by 
GMRES on step i is measured as |ri|, where fi is the residual, Vi = Axi — b. 
There holds that 

|fi| = inf \p{A)fo\, (5) 

pePi 
p(0) = l 

in which Pi is the space of monic polynomials of degree i. The matrix valued 
polynomial in the above minimization problem can be evaluated using Dunford 
integral mu- Let the open set V C C be such that cf{A) C U and dU is 
the union of rectifiable positively oriented Jordan curves. The set (j{A) is the 
spectrum of A. Application of the Dunford integral gives 

p{A) = -F( p{z){zl - A)~'^dz. (6) 

J'TTI Jqu 

This integral can be used to derive estimates for equation (|5|). Let A^ satisfy 
the assumptions made on the set U and in addition let 

Ae c Ae. 

This is, \zl — A\ < e~^ M z ^ dA^. In our case, Ag is an approximation for the 
pseudospectral set. Estimating the integral gives 

\p{A)ro\ < b(^)l|ro| < sup \p{z)\. (7) 

zeA. 

Gombining equations dU and 0 leads to the GMRES convergence estimate 



I^Ael 

27re 


inf sup \p{z)\. 
p^Pi zeK 

P(0) = 1 


(8) 


As we will illustrate in Section 4, this bound is useful for deriving worst case 
behavior of the GMRES convergence rate. 

The convergence bound ([8]) is meaningful only if one can solve the complex 
polynomial minimization problem. Typically, the set Ag is replaced with a 
larger set on which the minimization problem can be solved analytically. In 
simple cases, Ag can replaced with a circular or an elliptical domain, [2T]. Due 
to the constraint p(0) = 1, this approach gives useful information only when 
the circle or ellipsoid containing Ag does not contain the origin. When this is 
the case, one can try to apply so-called bratwurst shaped domains [15]. As 
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the name suggests, a bratwurst shaped domain can curl around the origin and 
it can be used to derive convergence estimates for the minimization problem. 
The bratwurst shaped domains can be applied with the inclusion and exclusion 
regions derived in this paper. However, the construction given in m is not 
simple, and cannot yield easy to use a priori bounds. 

In order to verify the analytically derived inclusion and exclusion results, 
we compute examples of the pseudospectral sets. Several different strategies for 
computing pseudospectrum have been proposed, see [ 23 ] and references therein. 
Several software packages, such as EigTool, are also freely available 0. 

To have full control over the computation of the pseudospectrum, we have 
chosen to use our own implementation of GRID - approach to compute pseu¬ 
dospectral sets. In the GRID-approach, a mesh is placed in the complex plane 
and the norm of the resolvent is computed for each grid point. The computed 
data is used to isolines describing the set. In the simplest case, the norm is 
computed as the largest singular value of the matrix {z — A)~^. Clearly, such 
an approach is very expensive for large number of points and large matrices. 
The process can be sped up by adapting the computational grid to the resolvent 
norm or by using a suitable matrix factorization to speed up the evaluation of 
the largest singular value. We have opted to speed up the computation by using 
an adaptive strategy to refine the computational grid. An initial triangular grid 
is placed in the complex plane. The grid is iteratively refined to conform to the 
shape of the resolvent norm. We use a refinement strategy based on splitting 
triangles intersecting with pre-specified level sets of the resolvent norm. This 
guarantees higher resolution at interesting regions of the complex plane. 

3 Abstract Framework 

In this section,we derive inclusion and exclusion regions for the pseudospectral 
set. For this purpose, it is easier to bound the complement of Ag, i.e. 

A^={|(^7-A)-i|<e-i}. (9) 

The inclusion and exclusion regions will lead to a set containing the pseudospec¬ 
trum. If the boundary of this set is a rectifiable Jordan curve, it can be used 
in connection with equation (jH|) to compute convergence estimates for the GM- 
RES method. The exclusion will be a disc around the origin. For the results to 
be meaningful, the exclusion should not be fully contained in the inclusion. If 
this is the case, the polynomial minimization problem in equation (|8|) does not 
tend to zero and the bound does not provide useful information. This has to be 
studied separately for each problem. 

When {zl — A) is non-singular, the matrix norm in equation (|9|) is defined 
as 

( 10 ) 

uEV \^u\ 

^see the Pseudospectral Gateway, http://www.cs.ox.ac.uk/pseudospectra/ 
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To eliminate the inverse and to establish a connection to the weak problem, we 
define an auxiliary vector Xy G such that 


(^zl A^Xy — 


( 11 ) 


Estimates for the resolvent norm are derived using the auxiliary variable. First, 
we establish the stability bound \xy\ < f{z)\xu\- When f{z) is bounded from 
above, this implies that {zl — A) is non-singular. In this case, the auxiliary 
vector is uniquely defined and 

Xy = {zl - A)~^Xu. 

The norm ([iQ]) can be estimated using the stability estimate for Xy as 

\(zI-A)-^\ = sup^§\<fiz). ( 12 ) 

uEV \^u\ 

We begin by taking advantage of the stability of the weak problem. Assumption 
12.11 Due to the duality between coefficient vectors and functions, stability of the 
weak problem implies stability of the linear system. As all finite dimensional 
norms are equal, there exists positive constants a, aw > 0 independent of u 
such that 

c^\xu\ < ll'^ll and aw\xu\ < ll'^llu^ G V. (13) 

When the derived framework is applied to a specific problem, a and aw are 
typically dependent on the mesh size. The dependency of these constants on 
relevant problem parameters are discussed in Section 4. Combining these norm 
equivalences with Assumption 12.11 leads to the following corollary. 

Corollary 3.1. Let b e and Xu he sueh that Axu = h. Then there holds that 

\Xu\<C2sK 

Where C 2 S •= • 

Proof. Let q e V he such that 

{q,v)w=xlb Vu G y. 


where (*, Ow is inner product on W. Using Cauchy-Schwarz inequality and the 
norm equivalence given in equation ([13]) there holds that \\q\\w < c^^\b\. Via 
this construction, vector b defines an antilinear functional on W' as L{v) := 
{q,v)w‘ By the definition of the dual norm and Cauchy-Schwarz inequality 


\\L\\w’ 


sup 

wew 


||w||vi/ 


< sup 
wew 


I|w||h/ 


IklIvK- 


(14) 


It follows that 

\\L\\w' < 

Combining the above equation with Assumption 12.11 and equation ([T3|) , we ob¬ 
tain 

aw(^\xu\ < Cs\b\. 


□ 
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The above Corollary essentially gives a lower bound for the smallest singular 
value of A. There holds that 

/ A\-l • 

crmin[A) = mm —^— 

\x\ 

so, that C 2 S ^ (^min{A). Corollarv 13.11 can be used to derive exclusion region 
near the origin. We give here a direct proof that fits well to the framework of 
the paper. Same result can be established from the lower bound for the smallest 
singular value by using Theorem 3 from m- 

Theorem 3.1. Let Assum.ntion A^. 1\ hold and let C 25 be as defined in Corollary 
ro Then there holds that 

BiO, -T-e)c K, 

^2S 

in whieh B{zo^r) :={ 2 :gC| \z — Zo\ <r }. 

Proof. From the definition of the auxiliary variable m it follows that 

AXy = ZXy — Xu 

Application of Corollary 13.11 gives 

\Xv\ < C 2 S {\z\\Xy\ + |f„|) 


i.e. 


\Xi; I ^ 


C 2 S 


1-C25k| 


(15) 


When | 2 )| < C 25 , the above bound implies that {zl — A) is non-singular. In this 
case, combining equations m and m gives 




To obtain the exclusion region, we set 

C 2 S 


1 - ^ 251^1 


< e 


which gives the bound 


U < 


C 2 S 


□ 


The inclusion is obtained by relating pseudospectrum to FOV. The following 
Theorem is proven e.g. in, m- For completeness, we give a proof using the 
notation used in this Section. 






Theorem 3.2. Let 5'^ := { z G C | dist{z,FOV{A)) <e} in which 


dist{z^ Q) := inf \z — q\. 
qeQ 


Then there holds that C 5'e. 

Proof. The auxiliary variable is defined as 

{A — Zl)Xy = 

Testing the above equation with Xy gives 


^AXy — ZX^Xy = XlXu- 

Using Cauchy-Schwarz inequality gives 

X%AXy 
XlXy 


> \xlAXy - ZXlXy \ = 


This is, 


A ^ 


ry* ^ rr* 


< \Xy 


By the definition of FOV(A) in equation (j4]) there holds that 

\xu\ > dist( 2 :,FOU(A))|f^|. 


□ 

Theorem 13.21 gives tools for deriving an inclusion for the pseudospectrum. 
The FOV is directly related to the boundedness properties of the sesquilinear 
form of the original problem. This relation arises from the connection x^Axy = 
a{v,v). The simplest estimate follows from boundedness of the sesquilinear 
form. Assume that there exists C > 0 such that 

\a{u,u)\<C\\u\\^ \/u&V. 

Then there holds that 

FOV {A) c 5(0, C). 

This is a very crude estimate, but it demonstrates how FOV can be bounded 
in simple cases. However, as we will see, more refined estimates are required to 
avoid inclusion of zero to the approximate pseudospectrum. 

4 Examples 

In this section, we demonstrate the presented theory with three examples. In all 
examples, we assume that UcM^,d = 2,3isa bounded domain with Lipschitz 
continuous boundary. We use standard notation for Sobolev spaces, see [1]. 
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The finite element space V is defined as 


V := {u e H^(n) I u e F^(K) VK e T } 


( 16 ) 


where T is a shape regular triangular or tetrahedral partition of [T]. This 
is, V is the space of first order Lagrange finite elements. The space of first 
order polynomials over set K is denoted by P^(K) and the mesh-size by h , 
respectively. 

The presented theoretical results are independent of the domain, but the 
actual numerical examples are computed on Cl = ( — 1,1)^ \ (0,1)^. The meshes 
used in the tests are generated from a coarse mesh with approximately 100 
nodes using uniform refinement. The coarse mesh is called level one mesh, once 
refined coarse mesh as a level two mesh and so on. 

Throughout this Section, c, C > 0 are generic positive constants independent 
of mesh size h, solution, load, and parameters of the weak problem, if not 
otherwise stated. They may depend on the shape regularity constant of the 
partition T and the domain Cl. 

4.1 Poisson equation 

We begin by considering the finite element discretization of the Poisson equation: 
Find G Vb such that 


(V^,Vu) = (/,u) Vug lb). 

In which Vq = V D Hq{CI) and / G L‘^{Cl). This is 

a{u,v) := (Vi4, Vu) and L{v) := (/, u) 
so that L G {L^{Cl)y. We use the standard iL^-norm 

ll'“lli := (Vu, Vw) + {u,u) 


(17) 


for the space Vo¬ 
lt is straightforward to see that the matrix A related to problem (pT|) is sym¬ 
metric and positive definite, [1]. The convergence of iterative methods for such 
linear systems can be analyzed using much easier techniques than pseudospec¬ 
trum. However, such a simple example is useful for demonstrating what kind 
of information on GMRES convergence can be obtained based on the inclusion 
and exclusion results. 

Pseudospectrum of a normal matrix can be easily computed from it’s eigen¬ 
values. All normal matrices are unitary diagonalizable, hence there exists a 
diagonal D G C^^^and a unitary Q G such that A = Q*DQ. Based on 

this expansion, there holds that 


\{zI-A) ^\ = \{z-D) = max K^-A) ^ 

' ' ' ' Aecr(A) 
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Thus, pseudospectrum of any normal matrix is a union of discs centered around 
it’s eigenvalues A^, 

V = utiS(Ai,e). 

The pseudospectrum for level one mesh is visualized in Figure ITT] for different 
values of e. 



Figure 1: Pseudospectral set for the Poisson problem on level one mesh. For 
sufficiently small e, the set is composed of disjoint disks with radius e. 

Next, we derive inclusion and exclusion regions using Theorem 13.11 and 13.21 
First, we need to establish a stability estimate satisfying Assumption 12 .11 As we 
are interested in mesh size explicit bounds, we use h-explicit norm equivalences 
instead of equation (p!3|) . For the Poisson problem, stability estimate follows 
from the weak problem (HZI by using Poincare-Friedrichs inequality. Let u eVq 
be the solution to dlZD then there exists a constant C > 0 such that 

Following this stability estimate, we choose the space W as and II • llvi/ = 

II • II 0 - Norm equivalences between L‘^{Q)- and the Euclidian norm can 

be derived in the finite element space V using the scaling argument and inverse 
inequality, m- There exists c and C such that 

ch^^‘^\xu\ < ||i^||o < Ch^^‘^\xu\ \/u eV (18) 

and 

ch^^‘^\xu\ < ll'i^lli < Ch^^‘^~^\xu\ G V. (19) 
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Now, we can apply Corollary 13.II to derive a stability constant for the linear 
system arising from the weak problem (ED- Let Xu be such that Axu = b. Then 
by Corollary 13.11 and the h-explicit norm equivalences, there exists a constant 
C such that 

\xu\ < Ch~‘^\b\. 

Application of Theorem 13.11 gives the following exclusion near the origin, 

B{0,Ch^-€) c A^. 

We proceed by deriving an inclusion for FOV (A), which together with Theorem 
13.21 gives inclusion for A^. It is easy to derive the estimates 

^a{u, u) = ^\\Vu\\l = 0 \/ueVo 


and 

ch'^lxuf < ^a{u,u) < Ch^~‘^\xu\‘^ G Vb- 
So that FOV {A) C { x G M | ch^ < x < Ch^~‘^ }. An application of Theorem 
[321 gives the inclusion C in which 

:= { 2: G C I dist (2;, { X G M I ch^ < x < Ch^~^ }) < e }. 


The above inclusion and exclusion regions give us an approximation of pseu¬ 
dospectrum Ag, 


Ae ■.= S,\B{0,Cih^ -e). 


Where the constant (7i > 0 is independent of h and e . To exclude the origin 
from this approximate pseudospectrum, we have to choose the parameter e as 
e < Cih^. In this case, the length of the boundary curve around the approximate 
pseudospectrum satisfies |9Ae| < C 2 h^~^ for some C 2 > 0 independent of h and 
e. When combined with equation ^ approximate pseudospectrum gives the 
GMRES convergence bound 


27re 


inf 

pcA 

p(0) = l 


sup \p{z)\\fo\ 

zeAe 


Ve < Cih'^ 


( 20 ) 


The set A^ can be covered either with an ellipsoid or a circle and the mini¬ 
mization problem can be solved using estimates given in [iniEi]. There holds 


that 


inf sup \p{z)\ < 

P^Pi zeB{c,r) 
p(0) = l 



i 


Although the estimate could be optimized with respect to parameter e, we have 
chosen e = 0.5(7ih^, which gives correct asymptotic behavior with respect to h. 
Using this e and c = ( 72 ^^“^, the circle based bound leads to the estimate 


Ini < 


C2h 


-2 


ttCi 


, ^ + 2 C ^2 . 


Ini 
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When the termination criteria for GMRES is chosen such that the relative resid¬ 
ual satisfies |ri||ro|“^ < to/, the above estimate gives the required number of 
iterations N as 

N « -^h-^ (^logtol - log ( 21 ) 

Our approximate pseudospectrum cannot capture the behavior of Ag for very 
small values of e. For example in the current case, the exact pseudospectrum is 
composed of small discs with boundary length 27re. Let e be such that the discs 
generating the pseudospectrum do not intersect. Any finite union of disjoint 
disks satisfies the conditions placed on the set U in the Dunford integral. Using 
equation (| 6 |) we obtain the estimate 

\p{A)\ <-^ F lp(^)l ^ ^0 sup \p{z)\, 

which is valid for sufficiently small e. Here Nq is the number of disjoint eigen¬ 
values of A. For quasi-uniform meshes, there exists C such that Nq < Ch~^ so 
that 

\p{A)\<Ch-‘^ sup|p(^)|. 

2:eAg 

Combining the above estimate with equation (j5j) gives 


< Ch-"^ inf 
ro| pePi 

p(0) = l 


sup \p{z)\. 

2;GAe 


( 22 ) 


This estimate based on the exact set Ag has a different multiplicative term 
in comparison to equation pn]). Interestinly, for d = 1, multiplicative term 
is smaller, for d = 2 it is equivalent and for d = 3 bigger. Regardless of 
the multiplicative constant, the estimate (| 22 | can deliver improved convergence 
number estimates. The best possible bound can be obtained at the limit e = 0, 
when the minimization problem can be solved using Chebychev polynomials, 
see e.g. [ 2 T]. Based on the FOV, the condition number n{A) < Ch~‘^. We 
obtain an estimate for the number of iterations 


N ^ —Ch ^(log(to/) — log(Cd ^)) (23) 

The main difference between the estimates m and (|23|l is in in the power 
of the mesh size h. For the particular problem, this difference is due to the fact, 
that the set Ag cannot capture the behaviour of the pseudospectrum for small e. 
For complicated problems, such knowledge is very difficult to come by and one 
has to be satisfied with worst case estimates, such as equation (|2T]). The second 
difference between the two estimates is in the additive terms. These additive 
terms are relevant only when tolerance is of the same order of magnitude with 
Ch~^^ which requires usage of very fine mesh sizes 
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4.2 Helmholtz equation with absorbing boundary condi¬ 
tions 

The Helmholtz equation with first-order absorbing boundary conditions is a 
more realistic example for the analysis presented in this paper. The weak prob¬ 
lem reads: Find u G such that 

a{u,v) = L{v) (24) 


in which 

a{u, v) := (Vu, Vv) + ifc (u, v)g^ - k^(u, v) and L{v) := (/, v) + {g, v)q^ . (25) 

The parameter k G k > 0, f G L^{Q) and g G L^{dCl). The inner product 
(*5 ')dn is the standard L^-inner product over dQ. The stability of this problem 
has been analyzed in domains excluding any resonant behavior, [19]. 

Theorem 4.1. Let Li be a bounded, star shaped domain with a smooth boundary 
and let u G be the solution to problem [2^ . Then there exists a eonstant 

Cs > 0 independent of u,f,g and n sueh that 

||w|U < (ll/llo + 115110 ,an), 

in whieh the norm || • \\^ is defined as 

\\u\\l:=\\Vu\\l + K^u\\l (26) 

The finite element approximation Uh is defined as: Find Uh ^ V such that 

a{uh, v) = (/, v) + {g, u)ao Vu G V. 

When the solution has 11)-regularity, the existence of a unique solution to 

this problem can be guaranteed, when the mesh size requirement n^h « 1 is 
satisfied, [T3[T1[T1]. In this case, there exists a constant C such that the a 
priori error estimate 


\\u - UhW^ < Ch (ll/llo + Hallo,an) • ( 27 ) 


holds. 

Due to the boundary term in {u,v)q^, problem (|M|) leads to a linear system 
with a non-normal coefficient matrix. As the boundary term depends on /^, it is 
complicated to determine if the non-normality is meaningful or not. In addition, 
due to the relation between the wave-number and the mesh size it is difficult to 
study the asymptotic behaviour of GMRES, when n tends to infinity. 

When the mesh size is sufficiently small so that the a priori error estimate 
holds. Theorem 14.11 implies stability of the discrete problem. We obtain 

ll'^/ilU ^ (1 + C'h)(||/||o + lbllo,o) < C'dI/llo + lbl|o,o)- (28) 

This discrete stability estimate holds under the following assumptions. 
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1 



Figure 2: The pseudospectral set for the Helmholtz equation with first order 
absorbing boundary conditions. The parameter k = Sn and level three mesh is 
used in the upper figure and level four in the lower one. 


Assumption 4.1. Assume that Q is a hounded, star shaped domain with a 
smooth boundary, the mesh size h is sueh that n?h « 1 and the solution u to 
problem [2^ has -regularity. 

Following the discrete stability result (|28|) we choose the space W as L^{Q) 
with the norm IHIw = IM|o- The h- explicit norm equivalences given in equation 
([T8|) can be used for this space. As we are interested in wavenumber and the 
mesh size explicit estimates, we use the /^-dependent norm given in equation 
([26|) for the space V. Norm equivalences for this /^-dependent norm are easily 
established using equation m and m as 

chvh^^‘^\xu\ < II^^IU ^ + Kh^^‘^)\xu\ \/u G V, (29) 

for some c, C. Application of Corollary 13.11 gives the stability estimate for the 
coefficient vector 

(30) 

rv 

Using Theorem 13.11 leads to the exclusion region 

H(0, Cnh^ - e) C A^ 
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Figure 3: Pseudospectral set for the Helmholtz equation with first order absorb¬ 
ing boundary conditions. The parameter k. = Sn. The mesh levels is three on 
left and four on right. One can observe the convergence of the set when mesh 
size tends to zero. 

around the origin. To obtain an inclusion, we again derive an inclusion for 
FOV{A) and apply Theorem I3.21 The sesquilinear form satisfies the bounded¬ 
ness estimates 

|5Ra(w,u)| < Cllwll^ VweV 

and 

0 < 'Aa{u^u) < CKh^~^\xuf \/u 

for some C. The estimate between the L^{dCl)- and Euclidian norm used in 
above is derived using identical techniques as used for proving inequality ([18]) . 
When Assumption 14.11 is satisfied, combining the two boundedness estimates 
leads to the inclusion 

FOV{A) c { z G C I |z| < C and 0 < } . 

This set contains the origin, so it cannot be used to derive GMRES convergence 
bounds. In this case, the presented theory is genuinely required to understand 
GMRES convergence. 

To validate the derived inclusion and exclusion regions, we have computed 
examples from the exact set for k = Sir using mesh levels three and four. 
The results are visualized in Eigures[2]and[3l Although, the L - shaped domain 
used in computations does not have smooth boundary nor -regularity, the 

actual pseudospectral set is in good agreement with our theoretical results. Most 
importantly, when e is sufficiently large, the pseudospectrum curls around the 
origin as predicted. Due to the solution having less that 11)-regularity, the 

requirement on the mesh size on L-shaped domain just takes the form « 1, 
for some a < 2, depending on regularity of the exact solution. 

The approximate pseudospectrum could also be used to to derive conver¬ 
gence estimate for GMRES method using Bratwurst shaped domains to solve 


16 














the minimization problem. However, as a preconditioner would always be ap¬ 
plied, the current case is not very interesting hence we do not proceed further 
with it. 

4.3 Shifted-Laplace preconditioned Helmholtz equation 

The analysis of inclusion and exclusion regions is more complicated, when a 
preconditioner is applied to speed up the convergence of the GMRES method. 
Several different preconditioners have been proposed for problem see e.g. 
[6] . We consider here the shifted-Laplace preconditioner [8] . This preconditioner 
is based on solving an auxiliary problem on each step of the iteration. The 
auxiliary problem is defined as: For a given u find Pu G V such that 

b{Pu,v)=xlxu VugF. (31) 

The sesquilinear form b in the above equation is given as 

b{u, v) = (Vi4, Vu) + m {u, v) + ia{u, u), 

in which k is as defined in Section 4.2 and a G M, cr > 0. This is, a loss term 
icr(i4, v) is added to the sesquilinear from defined in equation ([25]). The addition 
of the loss term leads to a stability estimate on the finite element space V 
independent of the mesh size. Choosing v = Pu in equation m gives 

b{Pu, Pu) = Xp^Xu Vu G V. 

Taking imaginary part leads to 

K\\Pu\\dQ + (T||-Plt||o = '^^PuXu- 

This is, 

cr||Pu||o < ^^PuXu- 

Now, using Cauchy-Schwarz inequality and norm equivalence ([T8|) gives 

\\Pu\\o<Ca-^\\u\\o yueV (32) 

for some C. The matrix form of the preconditioner is denoted as where 

Bij = b{(pj, (fi). Hence, the problem to be solved by the GMRES method is 

AB~^x = 5 , X = B~^x. 

The rationale behind using shifted-Laplace preconditioners is that when a 
sufficiently large loss term is added, the action of the preconditioner can be 
efficiently evaluated using a multigrid method, [7]. When applied directly to 
solve the original problem (|24|) . multigrid methods face two challenges, |4]. The 
standard smoothing iteration is not stable and the coarse grid correction has to 
be made on a sufficiently fine mesh. The introduction of a loss term has been 
analyzed in m for a problem with zero Dirichlet boundary conditions. In this 
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case, additional losses improve the multigrid solver by allowing the coarse grid 
correction to be made on a coarser mesh. The coarse grid depends on the loss 
term, hence there is a tradeoff between the number of GMRES iterations and 
the cost of applying the preconditioner. Typically, the loss parameter is chosen 
as cr = 0.5/i:^. For simplicity, we will consider here only the exact preconditioner. 
This gives good insight on what one can expect from the inexact case. 

As we will see, a shifted-Laplace preconditioner can eliminate the mesh size 
dependency from the pseudospectral set. This is, the inclusion and exclusion 
regions are independent of the applied mesh size. This is a desired property, 
as the mesh size dependency in the non-preconditioned case leads quickly to 
an unbearably large number of iterations. The exclusion regions will, however 
depend on the ratio of k and a. 

The shifted-Laplace preconditioner has been previously analyzed in [25] by 
estimating the location of the eigenvalues. The existing analysis is not explicit in 
a and does not take the non-normality into account. In addition, the previous 
work does not include the exclusion region around the origin, which we can 
obtain using Theorem 3.1. and the stability result given in equation (|28|). 

To study the shifted-Laplace preconditioner, we interpret the matrix AB~^ 
as the matrix form of the sesquilinear form a{Pu, u)^ where a{u, v) is as defined 
in equation (|2^ and P in equation (|3T]) . A suitable stability estimate for this 
sesquilinear form is established by the following Corollary. 

Corollary 4.1. Let u eV be sueh that 

a{Pu,v) = {f,v) 

In addition, let Assumption ^- 1\ be satisfied, 
independent ofu,f,K,h and a sueh that 

\xu\ < (^1 + 

Proof. Application of equation (|28)) gives 

II^^IU<^ll/llo. (34) 

It follows from definition (|3T]) that 

a{Pu, u) = — icr(Pi4, u). 

Combining above with equation (|33|) gives 

\xu\'^ = {f,u) A\a{Pu,u). (35) 

Using Cauchy-Schwarz inequality, estimate (|34|) and norm equivalence ([T8|) gives 

|^„|<C/i''/'(||/||o + C's^||/||o). 

□ 


Vi; G V. (33) 

Then there exists a constant C > 0 


-) 
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The above stability estimate is given in the norm ||i4|| = \xu\> Hence, we 
choose this as the norm of the space V. The above Corollary also suggest 
to choose the space W = L^{Q) as previously. With these choices, a direct 
application of Theorem 13.11 gives the exclusion 

H(0,C^-e)cA^ (36) 

When cr = 0, the preconditioner solves the problem exactly and = H(l,e). 
As the constant in above is C is independent of a and k, setting a = 0, leads 
to C < 1. A field of values based inclusion can be obtained as follows. There 
holds that 

= a{Pu,u) = xl^Xu — i(7{Pu,u). 

An inclusion for FOV follows by estimating the last term. By the stability result 
given in equation (|3^ and norm equivalence disi, there holds that 

(j{Pu,u) < (j\\Pu\\o\\u\\o < Cxl^Xu. 

This is, the FOV is located inside the set |1 — 2 :| < Ci. 

The polynomial minimization problem in the GMRES convergence bound 
® does not give any information on the convergence, when the approximate 
pseudospectrum is an annulus surrounding the origin. To apply the FOV based 
estimate, one has to explicitly know the constants in derived inclusion and 
exclusion regions to guarantee that this cannot happen. The constant Ci in 
the inclusion for FOV is related to the norm equivalence between L^{Q) and 
Euclidian norm. It is easy to see, that Ci = ^cond(M), where Mij = {(pi,(pj) 
is the mass matrix. In typical cases (7i ~ 4, so that derived inclusion is not 
useful when a = 0.5/^^ and the dimension of the exclusion tends to zero when n 
grows. 

Due to the close relation between the preconditioner and the original prob¬ 
lem, we can estimate the pseudospectrum using a problem specific technique. 

Lemma 4.1. There exists a positive eonstant C > 0 sueh that 


l^eC 


C 


( - - - 



C A 


c 

e 


Proof. There holds that A = A^ and B = B^. Using the identity \C\ = |C*| 
for any C G , it follows that 

\{zl - AB-^)-^x^\ \{zB - A)-^Bx^\ 

sup -p—j-= sup -p—- 

Xuec^ \^u\ xuec^ \^u\ 


Let Xy G be such that 


{A — zB)xv = Bxu‘ 
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As in Section 3, we establish the stability estimate \xv\ < f{z)\xu\> When f{z) 
is hnite, this estimate yields the desired bound. Testing with any x^ G 
gives 


(1 — z)a{v, w) — \(jz{v^ w) = ^Bxu- 

Assuming that 2 ; 7 ^ 1 and dividing by 1 — 2 ; yields 

, . iaz , ^ , u , icr , , 

a{v,w) - - - {v,w) = a(-- ,w) + -- {u,w). 

1 — z 1 — z 1 — z 

By adding an subtracting a suitable term, the above can be written as 

a(v - (1 - z)~'^u,w) - - (1 - z)~'^u,w) = . 2 ^u,w) 

1 — z (1 — Zf 


Choosing w — v — (1 — z) using the identity and taking 

imaginary part gives 

K\\v-{l-z)-^u\\l^Q^+a ^ ^y {u,v-{1-z)-^u) 

When 2 ; 7 ^ 1 and ^z — \z\^ > 0, this is 


1 

2 




the coefficient of the L?(yt) - term is positive and we obtain the estimate 

Using the norm equivalence given in equation ([18]) yields 
|f„ - (1 - z)-^Xu\ < ^ 

The stability estimate follows from the above equation and triangle inequality 
as 

IJ.I < lx. - (1 - x)-‘x.| + ^ < C (i^l + |x.|. 

□ 


To obtain an overview of the derived bounds we have computed the pseu¬ 
dospectrum for K = Idir and a = using the level three mesh. The 

results are presented in Figure IH Based on these results, analysis given in this 
Section seems to capture the behavior of the pseudospectrum rather well. In 
both cases, when e is sufficiently small, pseudospectrum is located inside the 
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disc , as predicted by Lemma [4.11 When the loss term is small, the 

pseudospectrum has a rather small diameter and is located close to 1. For 
large values of a, the set moves closer to the origin. These results are in good 
agreement with the exclusion given in equation (|36|) . 

The GMRES convergence bound gives usable information only if the origin is 
outside the approximate pseudospectrum. In the current case, this requirement 
limits the value of e and thus determines the GMRES convergence rate. We have 
studied the pseudospectrum close to the origin in more detail by using a bisection 
search to find x G M closest to the origin such that \{xl — = 2 • 10“^ 

for different hz between An and 647r for a = 0.5/i: and a = 0.5/^^. The results 
are visualized in Eigure [5l These results indicate, that the exclusion given in 
equation dssi corresponds well with the real behavior of the set. 




Eigure 4: Pseudospectrum for Example 4.3 with e = 1,10,100,1000. The pa¬ 
rameter K = IGtt and level three mesh was used. On left the loss term is chose 
as cr = 0.5/1: and on right as a = O.bK,^. The circle \) is visualized with a 
dashed line. 

A rigorous derivation of convergence estimate based on bratwurst shaped 
domains would require us to relate the parameters of these domains to A^, 
which is out of the scope of this paper. Our computations indicate that the 
pseudospectrum for sufficiently large e can be contained inside a circle, hence 
we will instead use the bound for circles to derive an approximate convergence 
rate. Based on the numerical and theoretical results, it seems to be reasonable 
to choose 

(J 

When (j = 0.5/i:^, there holds that Ag C 5(1,1— O.5/i:“^)05(O, e). To exclude the 
origin, we choose e = 0.25/i:“^. Using equation (jS]) and polynomial minimization 
over circles m, this leads to the estimate 

^ _ 

|ro| - Vl + 0.25«-i 
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Figure 5: The point x e R clos¬ 
est to the origin such that \{xl — 
A)~^\ = 2 • 10“^ for Example 4.3. 
The dependency is as predicted by 
the exclusion. Mesh level five was 
used in this computation. 


Figure 6: The number of GM- 
RES iterations required to solve the 
problem for different values of k. 
The stopping criteria was set as 
tol = 10“^. Level seven mesh was 
used in the computation. 


This, is the required number of iterations N to reach tolerance tol is 

N ^—Ahz log tol ^ Ahz log hz. (37) 

So, asymptotically, the dominating term is k. log hz. We cannot observe this effect 
in our numerical examples as it would require us to use extremely large values 
of K. For instance, when e = 10“^, k would need to be of the order 10^, before 
it has an impact on the required number of GMRES iterations. This means, 
that the non-normality is not practically relevant in this case. 

Estimate dSB rises the question, how should the stopping tolerance tol be 
chosen. Using the tools derived in this paper, the size of relative residual can 
be related to /^-dependent norm. As we have studied right preconditioning, the 
solution obtained from GMRES Xi = B~^Xi satisfies Vi = Axi — b = AB~^Xi — h. 
Hence, we can derive the estimate for the system without a preconditioner. The 
derived result holds for all left preconditioned systems. 

Lemma 4.2. Consider the problem Axu = b, in whieh A G and b G 

are related to the finite element diseretization of problem [2^ . Let Xu be sueh 
that \Axu — b\ <tol \b\. In addition, let Assumption \4A\ hold. Then there exists 
a eonstant C > 0 independent of tol,u,u, n and h sueh that 

\\u - u\\^ < Ctol (ll/llo + /i"^/^||5||o,£»n) . 

Proof. Denote r = Axu — b. There holds that Axu — b = A{xu — Xy)^ hence. 
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error e = Xu — Xu is a solution to the equation, 

Ae = r. 


Let the space W = L‘^{Q). Using the same construction as in the proof of 
Corollary 13.11 we define q gV such that {q,v)w = Vn G U and the linear 
functional L{w) = (g, w)w Using standard tools and the norm equivalence (p! 8 ]) 
gives ||L||w' < The stability estimate given in equation (|3Q|) leads 

to 

\\u-u\U<Ch-'^/M- 

Now, this can be written as 

11'^ — fi|U ^ Ch~^^^tol \b\. 


As there holds that 


\b\ 


b'^Xy 

max — 

\Xy\ 


{f,v) + {9,v)dn 

\Xv\ 


Cauchy-Schwarz inequality and norm equivivalence m gives 


M<Ch<^/^{\\f\\o + h-^/^g\\o) 


□ 

One should note that identical techniques that were used to prove the above 
Lemma can be used to derive a relation between the V - norm and tol for any 
finite dimensional variational problem satisfying Assumption 12.11 

We conclude by solving the shifted-Laplace preconditioned problem for right- 
hand side 

/ = exp (-103((a; - 0.5)2 ^ ^ 0 5 ) 2 ^^ 

and different values of k. The loss term for the preconditioner was chosen as 
a = and OAk and the level five mesh was used in the computations. The 

number of GMRES iterations is plotted in Figure [ 6 l In this case, we observe a 
linear relationship between k, and the number of iterations for a = 0.5k,^. The 
number of iterations stays constant for a = 0.5/^. These results are in good 
agreement with the estimate dSIl). 

5 Conclusions 

The main result of the paper is the derivation of exclusion region for pseu- 
dospectral set near the origin. Theorem 13.11 The derivation was made under 
Assumption 12.11 stability of the weak problem. All analysis was done a pri¬ 
ori, without constructing any matrices. Theorem 13.11 was applied in all three 
tests, and the derived results were in good agreement with the true behavior 
of the pseudospectral set. In addition, an inclusion region was derived using 
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the connection between FOV and the pseudospectrum. Boundedness estimates 
for FOV were derived based on the properties of the weak problem. All given 
analysis is applicable to a wide range of different problems. 

As demonstrated by the examples, the proposed inclusion and exclusion 
regions led to a worst case convergence estimate for the GMRES method. How¬ 
ever, the effect of this overestimation was significant only for extreme parameter 
values. As illustrated by the first example, more refined convergence estimate 
would require knowledge from behavior of pseudospectrum for e ^ 0. Such 
analysis is one direction for continuing this work. 

The aim of the paper was to investigate, if pseudospectrum based conver¬ 
gence estimate can be used for relating properties of weak form to convergence 
of GMRES. This was proven to be possible. As in Example 4.3, one needs to 
establish stability and boundedness of the preconditioned problem. The ap¬ 
plication of the derived theory will lead to inclusion and exclusion regions for 
pseudospectrum. Second possible direction for future work is to study differ¬ 
ent preconditioners and problems using the derived tools. Natural extension 
would be to investigate convergence of GMRES for time-harmonic Maxwell’s 
equations. 
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